library(data.table)
library(mefa)
library(doParallel)
library(emayili)
library(magrittr)
library(ggplot2)
library(nleqslv)


source('helper_functions.r')



runCounterfactualsEquity <- function(estimation) {
  
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  baseline <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  SCALE = 1710 / baseline[variable == 'total.volume']$baseline
  
  # Conforming loan limits
  cr.30 <- runCounterfactual(estimation = estimation,CR = 0.030,LIM = NULL,RGSE = NULL,NAME = '3pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE =0.25)
  cr.45 <- runCounterfactual(estimation = estimation,CR = 0.045,LIM = NULL,RGSE = NULL,NAME = '45pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE = 0.25)
  cr.60 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE = 0.25)
  cr.75 <- runCounterfactual(estimation = estimation,CR = 0.075,LIM = NULL,RGSE = NULL,NAME = '75pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE = 0.25)
  cr.90 <- runCounterfactual(estimation = estimation,CR = 0.090,LIM = NULL,RGSE = NULL,NAME = '9pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE = 0.25)
  cr.12 <- runCounterfactual(estimation = estimation,CR = 0.120,LIM = NULL,RGSE = NULL,NAME = '12pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,EQUITY_ISSUANCE_RATE = 0.25)
  cr = merge(cr.30,merge(cr.45,merge(cr.60,merge(cr.75,merge(cr.90,cr.12,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  # define rel columns
  rel.columns = c('rate.conforming','rate.jumbo','rate.spread','lender.profit','bank.profit','sb.profit','overall.surplus','individual.surplus','top.market.surplus','bottom.market.surplus','individual.surplus.top','individual.surplus.low')
  cr[variable %in% rel.columns,`3pct` := `3pct` - baseline ]
  cr[variable %in% rel.columns,`45pct` := `45pct` - baseline ]
  cr[variable %in% rel.columns,`75pct` := `75pct` - baseline ]
  cr[variable %in% rel.columns,`9pct` := `9pct` - baseline ]
  cr[variable %in% rel.columns,`12pct` := `12pct` - baseline ]
  cr[variable %in% rel.columns, baseline := 0 ]
  
  
  save(list = ls(),file = paste0('counterfactual_results/',estimation,'_counterfactuals_entry_no_exit_equity_jan_2021.rsave'))  
}

runCounterfactualsJumboSec <- function(estimation) {
  
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  baseline <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  SCALE = 1710 / baseline[variable == 'total.volume']$baseline
  
  # Conforming loan limits
  cr.30 <- runCounterfactual(estimation = estimation,CR = 0.030,LIM = NULL,RGSE = NULL,NAME = '3pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,JUMBO_SEC_RATE = 0.25)
  cr.45 <- runCounterfactual(estimation = estimation,CR = 0.045,LIM = NULL,RGSE = NULL,NAME = '45pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,JUMBO_SEC_RATE = 0.25)
  cr.60 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)#,JUMBO_SEC_RATE = 0.25)
  cr.75 <- runCounterfactual(estimation = estimation,CR = 0.075,LIM = NULL,RGSE = NULL,NAME = '75pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,JUMBO_SEC_RATE = 0.25)
  cr.90 <- runCounterfactual(estimation = estimation,CR = 0.090,LIM = NULL,RGSE = NULL,NAME = '9pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,JUMBO_SEC_RATE = 0.25)
  cr.12 <- runCounterfactual(estimation = estimation,CR = 0.120,LIM = NULL,RGSE = NULL,NAME = '12pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,JUMBO_SEC_RATE = 0.25)
  cr = merge(cr.30,merge(cr.45,merge(cr.60,merge(cr.75,merge(cr.90,cr.12,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  # define rel columns
  rel.columns = c('rate.conforming','rate.jumbo','rate.spread','lender.profit','bank.profit','sb.profit','overall.surplus','individual.surplus','top.market.surplus','bottom.market.surplus','individual.surplus.top','individual.surplus.low')
  cr[variable %in% rel.columns,`3pct` := `3pct` - baseline ]
  cr[variable %in% rel.columns,`45pct` := `45pct` - baseline ]
  cr[variable %in% rel.columns,`75pct` := `75pct` - baseline ]
  cr[variable %in% rel.columns,`9pct` := `9pct` - baseline ]
  cr[variable %in% rel.columns,`12pct` := `12pct` - baseline ]
  cr[variable %in% rel.columns, baseline := 0 ]
  
 
  save(list = ls(),file = paste0('counterfactual_results/',estimation,'_counterfactuals_entry_no_exit_jumbo_sec_jan_2021.rsave'))  
}


runCounterfactualsMain <- function(estimation) {
  
  
  
  print('doing it!')
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  baseline <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  SCALE = 1710 / baseline[variable == 'total.volume']$baseline
  
  # Conforming loan limits
  cr.30 <- runCounterfactual(estimation = estimation,CR = 0.030,LIM = NULL,RGSE = NULL,NAME = '3pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr.45 <- runCounterfactual(estimation = estimation,CR = 0.045,LIM = NULL,RGSE = NULL,NAME = '45pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr.60 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr.75 <- runCounterfactual(estimation = estimation,CR = 0.075,LIM = NULL,RGSE = NULL,NAME = '75pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr.90 <- runCounterfactual(estimation = estimation,CR = 0.090,LIM = NULL,RGSE = NULL,NAME = '9pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr.12 <- runCounterfactual(estimation = estimation,CR = 0.120,LIM = NULL,RGSE = NULL,NAME = '12pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  cr = merge(cr.30,merge(cr.45,merge(cr.60,merge(cr.75,merge(cr.90,cr.12,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  # define rel columns
  rel.columns = c('rate.conforming','rate.jumbo','rate.spread','lender.profit','bank.profit','sb.profit','overall.surplus','individual.surplus','top.market.surplus','bottom.market.surplus','individual.surplus.top','individual.surplus.low')
  cr[variable %in% rel.columns,`3pct` := `3pct` - baseline ]
  cr[variable %in% rel.columns,`45pct` := `45pct` - baseline ]
  cr[variable %in% rel.columns,`75pct` := `75pct` - baseline ]
  cr[variable %in% rel.columns,`9pct` := `9pct` - baseline ]
  cr[variable %in% rel.columns,`12pct` := `12pct` - baseline ]
  cr[variable %in% rel.columns, baseline := 0 ]
  
  # QE
  qe.m.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -1.00,NAME = 'm100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.m.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.25,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.m.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.10,NAME = 'm10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.base  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.p.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.10,NAME = 'p10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.p.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.25,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe.p.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  1.00,NAME = 'p100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  qe = merge(qe.m.100,merge(qe.m.25,merge(qe.m.10,merge(qe.base,merge(qe.p.10,merge(qe.p.25,qe.p.100,by='variable'),by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  qe[variable %in% rel.columns, m100 := m100 - baseline ]
  qe[variable %in% rel.columns, m25 := m25 - baseline ]
  qe[variable %in% rel.columns, m10 := m10 - baseline ]
  qe[variable %in% rel.columns, p10 := p10 - baseline]
  qe[variable %in% rel.columns, p25 := p25 - baseline ]
  qe[variable %in% rel.columns, p100 := p100 - baseline ]
  qe[variable %in% rel.columns, baseline := 0 ]
  
  
  # Limits
  lim.m.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 0.75,RGSE = NULL,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim.base   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim.p.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 1.25,RGSE = NULL,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim.p.none <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 5   ,RGSE = NULL,NAME = 'none',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim.lower  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_low',RGSE = NULL,NAME = 'all_low',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim.upper  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_high',RGSE = NULL,NAME = 'all_high',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  lim = merge(lim.m.25,merge(lim.base,merge(lim.p.25,merge(lim.p.none,merge(lim.lower,lim.upper,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  lim[variable %in% rel.columns, m25 := m25 - baseline ]
  lim[variable %in% rel.columns, p25 := p25 - baseline ]
  lim[variable %in% rel.columns, all_low := all_low - baseline ]
  lim[variable %in% rel.columns, none := none - baseline ]
  lim[variable %in% rel.columns, all_high := all_high - baseline]
  lim[variable %in% rel.columns, baseline := 0 ]
  
  
  save(list = ls(),file = paste0('counterfactual_results/',estimation,'_counterfactuals_entry_no_exit_jan_2021.rsave'))  
}



runCounterfactualsMain <- function(estimation) {
  
  
  
  print('doing it!')
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  b.2010 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2010)
  b.2011 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2011)
  b.2012 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2012)
  b.2013 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2013)
  b.2014 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2014)
  b.2015 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2015)
  b.2016 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2016)
  b.2017 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,CF_YEARS = 2017)
  
  
}


runCounterfactualsMain_by_Risk <- function(estimation,risk = 'LOW') {
  
  
  
  print('doing it!')
  
  ENTRY = T
  BANK_EXIT = F
  
  # Get baseline scale
  baseline <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT)
  SCALE = 1710 / baseline[variable == 'total.volume']$baseline
  
  # Conforming loan limits
  cr.30 <- runCounterfactual(estimation = estimation,CR = 0.030,LIM = NULL,RGSE = NULL,NAME = '3pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr.45 <- runCounterfactual(estimation = estimation,CR = 0.045,LIM = NULL,RGSE = NULL,NAME = '45pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr.60 <- runCounterfactual(estimation = estimation,CR = NULL ,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr.75 <- runCounterfactual(estimation = estimation,CR = 0.075,LIM = NULL,RGSE = NULL,NAME = '75pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr.90 <- runCounterfactual(estimation = estimation,CR = 0.090,LIM = NULL,RGSE = NULL,NAME = '9pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr.12 <- runCounterfactual(estimation = estimation,CR = 0.120,LIM = NULL,RGSE = NULL,NAME = '12pct',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  cr = merge(cr.30,merge(cr.45,merge(cr.60,merge(cr.75,merge(cr.90,cr.12,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  # define rel columns
  rel.columns = c('rate.conforming','rate.jumbo','rate.spread','lender.profit','bank.profit','sb.profit','overall.surplus','individual.surplus','top.market.surplus','bottom.market.surplus','individual.surplus.top','individual.surplus.low')
  cr[variable %in% rel.columns,`3pct` := `3pct` - baseline ]
  cr[variable %in% rel.columns,`45pct` := `45pct` - baseline ]
  cr[variable %in% rel.columns,`75pct` := `75pct` - baseline ]
  cr[variable %in% rel.columns,`9pct` := `9pct` - baseline ]
  cr[variable %in% rel.columns,`12pct` := `12pct` - baseline ]
  cr[variable %in% rel.columns, baseline := 0 ]
  
  # QE
  qe.m.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -1.00,NAME = 'm100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.m.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.25,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.m.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = -0.10,NAME = 'm10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.base  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.p.10  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.10,NAME = 'p10',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.p.25  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  0.25,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe.p.100 <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE =  1.00,NAME = 'p100',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  qe = merge(qe.m.100,merge(qe.m.25,merge(qe.m.10,merge(qe.base,merge(qe.p.10,merge(qe.p.25,qe.p.100,by='variable'),by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  qe[variable %in% rel.columns, m100 := m100 - baseline ]
  qe[variable %in% rel.columns, m25 := m25 - baseline ]
  qe[variable %in% rel.columns, m10 := m10 - baseline ]
  qe[variable %in% rel.columns, p10 := p10 - baseline]
  qe[variable %in% rel.columns, p25 := p25 - baseline ]
  qe[variable %in% rel.columns, p100 := p100 - baseline ]
  qe[variable %in% rel.columns, baseline := 0 ]
  
  
  # Limits
  lim.m.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 0.75,RGSE = NULL,NAME = 'm25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim.base   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = NULL,RGSE = NULL,NAME = 'baseline',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim.p.25   <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 1.25,RGSE = NULL,NAME = 'p25',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim.p.none <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 5   ,RGSE = NULL,NAME = 'none',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim.lower  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_low',RGSE = NULL,NAME = 'all_low',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim.upper  <- runCounterfactual(estimation = estimation,CR = NULL,LIM = 'all_high',RGSE = NULL,NAME = 'all_high',scale = SCALE,entry = ENTRY,ALLOW_BANK_EXIT = BANK_EXIT,risk = risk)
  lim = merge(lim.m.25,merge(lim.base,merge(lim.p.25,merge(lim.p.none,merge(lim.lower,lim.upper,by='variable'),by='variable'),by='variable'),by='variable'),by='variable')
  
  lim[variable %in% rel.columns, m25 := m25 - baseline ]
  lim[variable %in% rel.columns, p25 := p25 - baseline ]
  lim[variable %in% rel.columns, all_low := all_low - baseline ]
  lim[variable %in% rel.columns, none := none - baseline ]
  lim[variable %in% rel.columns, all_high := all_high - baseline]
  lim[variable %in% rel.columns, baseline := 0 ]
  
  
  save(list = c('cr','qe','lim'),file = paste0('counterfactual_results/',estimation,'_',risk,'_counterfactuals_entry_no_exit_jan_2021.rsave'))  
}


runCounterfactual <- function(estimation = 'estimation_16',CR = NULL,LIM = NULL,RGSE = NULL,EQUITY_ISSUANCE_RATE = 999,JUMBO_SEC_RATE = 999, NAME = 'baseline',scale = NULL,entry = F,ALLOW_BANK_EXIT = F,risk = 'ALL',CF_YEARS = 2017) {
  cf.data <- prepCounterfactualData(estimation,CR,LIM,RGSE,EQUITY_ISSUANCE_RATE = EQUITY_ISSUANCE_RATE,JUMBO_SEC_RATE = JUMBO_SEC_RATE,ALLOW_BANK_EXIT = ALLOW_BANK_EXIT)
  if(risk == 'LOW') {
    cf.data$market.data <- cf.data$market.data[FICO == 'l']
  } else if(risk == 'HIGH') {
    cf.data$market.data <- cf.data$market.data[FICO == 'h']
  }
  
  markets.data <- cf.data$market.data[YEAR %in% CF_YEARS]
  markets.data[,max.xi := max(abs(xi)),by='MARKET_ID']
  markets.data <- markets.data[max.xi < 10 ]
  markets.data[,MARKET_INCOME_BUCKET := 'middle']
  markets.data[market.log_income.mean < quantile(markets.data$market.log_income.mean,.25),MARKET_INCOME_BUCKET := 'low']
  markets.data[market.log_income.mean > quantile(markets.data$market.log_income.mean,.75),MARKET_INCOME_BUCKET := 'high']
  
  if(entry) {
    solved <- getAllEquilbriumInterestRatesWithEntry(markets.data = markets.data,linear.model = cf.data$linear.model,raw.draws = cf.data$raw.draws,demand.parameters = cf.data$demand.pars)  
  } else {
    solved <- getAllEquilbriumInterestRates(markets.data = markets.data,linear.model = cf.data$linear.model,raw.draws = cf.data$raw.draws,demand.parameters = cf.data$demand.pars)  
  }
  
  
  # Merge it with info about lenders
  ll <- merge(solved,markets.data[,c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','j','LENDER','JUMBO','HELD','MARKET_SIZE','MARKET_SHARE_DATA','MARKET_INCOME_BUCKET','ALTERNATIVE')],by=c('MARKET_ID','j'))
  
  # Get total size for adjustment
  if(!is.null(scale)) {
    ll[,SHARE_MODEL := SHARE_MODEL * scale]
  }


  
  # Stats
  volumes.m = ll[,j=list(total.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE,na.rm=T)/1e9,
                       conforming.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T)/1e9,
                       jumbo.volume = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'j'),na.rm=T)/1e9,
                       bank.volume  = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (LENDER == 'b'),na.rm=T)/1e9)]
  
  financing.m  = ll[,j=list(balance.sheet.financing = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * HELD,na.rm=T)/1e9,
                            balance.sheet.financing.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * HELD,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE,na.rm=T),
                            conforming.balance.sheet.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c') * HELD,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T),
                            shadow.bank.market.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (LENDER != 'b') ,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE ,na.rm=T),
                            shadow.bank.conforming.share = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c') * (LENDER != 'b') ,na.rm=T) / sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (JUMBO == 'c'),na.rm=T),
                            alternative_financing = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * ALTERNATIVE,na.rm=T)/1e9)]
  
  rates.m = ll[,j=list(rate.conforming = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'c'),na.rm=T),
                       rate.jumbo = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'j'),na.rm=T),
                       rate.spread = .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'j'),na.rm=T) -  .01 * weighted.mean(RATE_INST,w = MARKET_SIZE * SHARE_MODEL * (JUMBO == 'c'),na.rm=T))]
  
  profits.m = ll[,j=list(lender.profit = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)),na.rm=T)/1e9,
                         bank.profit   = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)) * (LENDER == 'b')/1e9,na.rm=T),
                         sb.profit   = sum(LOAN_SIZE * SHARE_MODEL * MARKET_SIZE * (RATE_INST - MC)/100 * ( (1- (1+0.05)^(-10) ) / (0.05)) * (LENDER != 'b')/1e9,na.rm=T),
                         mean.nfsb   = weighted.mean(N_LENDERS,w=as.numeric(LENDER == 'n'),na.rm=T))]
  
  # Surplus is in units of surplus per dollar of loan...
  # First, drop the repeats (should be equal by market)
  ll.surplus <- ll[!duplicated(MARKET_ID)]
  surplus.m = ll.surplus[,j=list(overall.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * LOAN_SIZE,na.rm=T)/1e9,
                         individual.surplus = weighted.mean(SURPLUS * SHARE_MODEL * LOAN_SIZE,w=,na.rm = T),
                         top.market.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * (MARKET_INCOME_BUCKET == 'high') * LOAN_SIZE,na.rm=T )/1e9,
                         bottom.market.surplus = sum(SURPLUS * SHARE_MODEL * MARKET_SIZE * (MARKET_INCOME_BUCKET == 'low') * LOAN_SIZE,na.rm=T )/1e9,
                         individual.surplus.top = weighted.mean(SURPLUS_HIGH * SHARE_MODEL * LOAN_SIZE,na.rm=T),
                         individual.surplus.low = weighted.mean(SURPLUS_LOW * SHARE_MODEL * LOAN_SIZE,na.rm=T))]
  
  all.columns = cbind(volumes.m,financing.m,rates.m,profits.m,surplus.m)
  result <- melt(all.columns)
  names(result) = c('variable',NAME)
  
  return(result)
}

getAllEquilbriumInterestRates <- function(markets.data,linear.model,raw.draws,demand.parameters) {
  
  ## Parallelized and non-parallelized
  all.markets <- unique(markets.data$MARKET_ID)
  
  CORES = 30
  
  ## Parallelized
  if(CORES > 1) {
    registerDoParallel(cores = CORES)
    solved.model <- foreach(ii = 1:length(all.markets), .combine = rbind) %dopar% {
        current.market.id <- all.markets[ii]
        current.market.data       <- markets.data[MARKET_ID == current.market.id]
        getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters)
    }
  } else {  
    for(ii in 1:length(all.markets)) {
      current.market.id <- all.markets[ii]
      print(current.market.id)
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      if(ii == 1) {
        solved.model <- getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters)
      } else {
        solved.model <- rbind(solved.model,getEquilibriumInterestRate(current.market.data,linear.model,raw.draws,demand.parameters))
      }
      
    } 
  }
  
  
  return(solved.model) 
  
}


getAllEquilbriumInterestRatesWithEntry <- function(markets.data,linear.model,raw.draws,demand.parameters) {
  
  ## Parallelized and non-parallelized
  all.markets <- unique(markets.data$MARKET_ID)
  
  CORES = 30
  
  ## Parallelized
  if(CORES > 1) {
    registerDoParallel(cores = CORES)
    solved.model <- foreach(ii = 1:length(all.markets), .combine = rbind) %dopar% {
      current.market.id <- all.markets[ii]
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters)
    }
  } else {  
    for(ii in 1:length(all.markets)) {
      current.market.id <- all.markets[ii]
      print(current.market.id)
      current.market.data       <- markets.data[MARKET_ID == current.market.id]
      if(ii == 1) {
        solved.model <- getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters)
      } else {
        solved.model <- rbind(solved.model,getEquilibriumInterestRateWithEntry(current.market.data,linear.model,raw.draws,demand.parameters))
      }
      
    } 
  }
  
  
  return(solved.model) 
  
}

getEquilibriumInterestRate <- function(market.data,linear.model,raw.draws,demand.parameters) {
  
  market.data <- market.data[order(j)]
  transformed.draws <- transformRawDraws(raw.draws,market.data,demand.parameters)
  
  objective.function <- function(x,info = F) {
    
    # Assume you get 'x' in the right order
    market.data$RATE_INST <- x
    
    # Get markup
    mr <- calculateMarginalRevenuesAndCounterfactualShares(market.data = market.data,transformed.draws = transformed.draws,pars = demand.parameters,linear.model = linear.model)
    mr <- merge(mr,market.data[,c('j','RATE_INST','MC','N_LENDERS')],by='j')
    mr[,foc := MC + MARKUP_MODEL - RATE_INST] 
    val <- mr$foc
    if(info) {
      return(mr[,c('MARKET_ID','j','SHARE_MODEL','LOAN_SIZE','MARKUP_MODEL','RATE_INST','MC','foc','SURPLUS','SURPLUS_HIGH','SURPLUS_LOW','N_LENDERS'),with=F])
    } else {
      return(val)
    }
  }
  
  
  soln <- nleqslv(x = market.data$RATE_INST,fn = objective.function)
  result <- objective.function(soln$x,info = T)  
  return(result)  
}


getEquilibriumInterestRateWithEntry <- function(market.data,linear.model,raw.draws,demand.parameters) {
  
  market.data <- market.data[order(j)]
  transformed.draws <- transformRawDraws(raw.draws,market.data,demand.parameters)
  
  cost_m = 10.37
  cost_s = 7.30
  cost_N = 815
  
  
  objective.function <- function(x,info = F) {
    
    # Assume you get 'x' in the right order
    toWork <- market.data
    
    
    toWork$RATE_INST <- x[1:nrow(market.data)]
    toWork[LENDER == 'n',N_LENDERS := x[1 + nrow(market.data)]]
    jj = toWork[LENDER == 'n']$j
    
    # Get markup
    mr <- calculateMarginalRevenuesAndCounterfactualShares(market.data = toWork,transformed.draws = transformed.draws,pars = demand.parameters,linear.model = linear.model)
    mr <- merge(mr,toWork[,c('j','RATE_INST','MC','MARKET_SIZE','N_LENDERS')],by='j')
    mr[,foc := MC + MARKUP_MODEL - RATE_INST] 
    mr[,VARIABLE_PROFIT := MARKET_SIZE * SHARE_MODEL/N_LENDERS * LOAN_SIZE * (MARKUP_MODEL / 100)/1000 ] # in k dollars, per lender]
    nfsb.var.profit = max(0,sum(mr[j %in% jj]$VARIABLE_PROFIT))
    nfsb.E_LENDERS = cost_N * pnorm(log(nfsb.var.profit),cost_m,cost_s)
    D_LENDERS =  nfsb.E_LENDERS - x[1 + nrow(market.data)]
    val <- c(mr$foc,D_LENDERS)
    if(info) {
      return(mr[,c('MARKET_ID','j','SHARE_MODEL','LOAN_SIZE','MARKUP_MODEL','RATE_INST','MC','foc','SURPLUS','SURPLUS_HIGH','SURPLUS_LOW','N_LENDERS'),with=F])
    } else {
      return(val)
    }
  }
  
  
  soln <- nleqslv(x = c(market.data$RATE_INST,market.data[LENDER == 'n']$N_LENDERS[1]),fn = objective.function)
  result <- objective.function(soln$x,info = T)  
  return(result)  
}



prepCounterfactualData <- function(estimation = 'estimation_5',CR = NULL,LIM = NULL,RGSE = NULL,JUMBO_SEC_RATE = 999, EQUITY_ISSUANCE_RATE = 999,ALLOW_BANK_EXIT = F,
                                   override_A = NULL,override_B = NULL, override_S = NULL,override_linear_demand = NULL,
                                   override_linear_year = NULL,override_linear_supply = NULL, override_nonlinear_supply = NULL) {
  
  
  
  # Load demand estimation data
  demand.estimation <- loadEstimation(estimation,override_A = override_A,override_B = override_B,override_S = override_S, override_linear_demand = override_linear_demand)
  market.data       <- demand.estimation$market.data
  linear.model      <- demand.estimation$linear.model
  demand.pars       <- demand.estimation$pars
  raw.draws         <- demand.estimation$raw.draws
  
  
  # Counterfactual conforming loan limit?
  if(!is.null(LIM)) {
    if(LIM == 'all_low') {
      market.data[, market.conforming_loan_limit := 417000]
    } else if(LIM == 'all_high') {
      market.data[, market.conforming_loan_limit := 625000]
    } else {
      market.data[, market.conforming_loan_limit := market.conforming_loan_limit * LIM]
    }
  }
  
  # Load supply estimation data
  supply.estimation <- loadSupplyEstimationForCF(estimation,override_linear_year = override_linear_year,
                                                 override_linear_supply = override_linear_supply,
                                                 override_nonlinear_supply = override_nonlinear_supply)
  estimation.data <- loadSupplyEstimationData(estimation,demand.estimation)
  
  # Counterfactual GSE costs?
  if(!is.null(RGSE)) { 
    supply.estimation[,MC_NONLIN_GSE := MC_NONLIN_GSE + RGSE]
  }
  
  ECAP_THRESHHOLD = 0.01
  
  # Get marginal costs for banks
  bank.data       <- fread('../data/final/capitalization_data_sep9.csv')
  merged.bank <- merge(bank.data,supply.estimation[LENDER == 'b',c('YEAR','LENDER','JUMBO','FICO','PURPOSE','MC_LIN_YEAR','MC_LIN_RISK_JUMBO','MC_LIN_LABOR','MC_NONLIN_GSE','MC_NONLIN_RE','MC_NONLIN_PSI'),with=F],by=c('YEAR','JUMBO','PURPOSE'),allow.cartesian = T)
  #merged.bank[,ACTIVE_BASELINE := (BANK_CR - rho_hat) > ECAP_THRESHHOLD]
  merged.bank[,ACTIVE_BASELINE := 1]
  
  if(!is.null(CR)) {
    merged.bank[,rho_hat := CR]
  }
  
  if(ALLOW_BANK_EXIT) {
    merged.bank[,ACTIVE_CF := (BANK_CR - rho_hat) > ECAP_THRESHHOLD]
  } else {
    merged.bank[,ACTIVE_CF := 1]
  }
  

  
  
  merged.bank[,ECAP_BANK := pmax(ECAP_THRESHHOLD,BANK_CR - rho_hat)]
  merged.bank[,mc_jumbo_nonlin_balance      := pmin(BANK_CR^2 * MC_NONLIN_RE * MC_NONLIN_PSI * xi_j * (ECAP_BANK)^(-1 * (MC_NONLIN_PSI + 1)),20)] # Just don't let anyone go crazy---they'd drop out anyway...
  merged.bank[,mc_conforming_nonlin_balance := BANK_CR^2 * MC_NONLIN_RE * MC_NONLIN_PSI * xi_g * (ECAP_BANK)^(-1 * (MC_NONLIN_PSI + 1))]
  merged.bank[,mc_conforming_nonlin_gse          := MC_NONLIN_GSE]
  merged.bank[JUMBO == 'j',mc_nonlin := pmin(mc_jumbo_nonlin_balance,mc_conforming_nonlin_gse + JUMBO_SEC_RATE,mc_conforming_nonlin_gse + EQUITY_ISSUANCE_RATE)] 
  merged.bank[JUMBO == 'c',mc_nonlin := pmin(mc_conforming_nonlin_balance,mc_conforming_nonlin_gse)]
  merged.bank[JUMBO == 'j',held := 1] # If you issue equity you're still holding it on balance sheet.....
  merged.bank[JUMBO == 'j',alternative := 0]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + JUMBO_SEC_RATE) < mc_jumbo_nonlin_balance,alternative := 1]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + EQUITY_ISSUANCE_RATE) < mc_jumbo_nonlin_balance,alternative := 1]
  merged.bank[JUMBO == 'j' & (mc_conforming_nonlin_gse + JUMBO_SEC_RATE) < mc_jumbo_nonlin_balance,held := 0] # Jumbo securitization
  merged.bank[JUMBO == 'c',held := 0]
  merged.bank[JUMBO == 'c' & mc_conforming_nonlin_balance < mc_conforming_nonlin_gse,held := 1]
  merged.bank[,MC := MC_LIN_YEAR + MC_LIN_RISK_JUMBO + MC_LIN_LABOR + mc_nonlin]
  bank.mc <- merged.bank[,j=list(LENDER = 'b',PCT_REMAIN = weighted.mean(ACTIVE_CF,w=N_ORIG) / weighted.mean(ACTIVE_BASELINE,w=N_ORIG),
                                 HELD = weighted.mean(held,w=N_ORIG * ACTIVE_CF),
                                 ALTERNATIVE = weighted.mean(alternative,w=N_ORIG * ACTIVE_CF),
                                 MC   = weighted.mean(MC,  w=N_ORIG * ACTIVE_CF)),by=c('YEAR','CBSA','PURPOSE','FICO','JUMBO')]
  
  # Combine it with fintech and non-fintech
  #non-bank
  nb.supply <- supply.estimation[LENDER != 'b']
  nb.supply[,HELD := 0]
  nb.supply[,PCT_REMAIN := 1]
  nb.supply[,ALTERNATIVE := 0]
  nb.supply[,MC := MC_LIN_LABOR + MC_LIN_YEAR + MC_NONLIN_GSE + MC_LIN_RISK_JUMBO]
  
  nb <- merge(estimation.data[LENDER != 'b',c('YEAR','CBSA','PURPOSE','FICO','JUMBO','LENDER')],nb.supply[,c('YEAR','FICO','PURPOSE','LENDER','HELD','ALTERNATIVE','PCT_REMAIN','MC')],by=c('YEAR','FICO','LENDER','PURPOSE'),allow.cartesian = T)
  

  combined <- rbind(nb,bank.mc)
  combined <- combined[order(YEAR,CBSA,PURPOSE,FICO,LENDER,JUMBO)]
  combined[is.na(MC),MC := 100]
  combined[is.na(PCT_REMAIN), PCT_REMAIN := 0]
  combined[PCT_REMAIN == 0,PCT_REMAIN := .01]
  
  capitalization.data <- combined
  
  supply.data <- merge(capitalization.data,market.data[,c('YEAR','CBSA','PURPOSE','FICO','MARKET_ID','JUMBO','LENDER','j','N_LENDERS')],by=c('YEAR','CBSA','PURPOSE','FICO','JUMBO','LENDER'))
  supply.data[,N_LENDERS := N_LENDERS * PCT_REMAIN]
  supply.data <- supply.data[,c('MARKET_ID','YEAR','CBSA','FICO','PURPOSE','j','LENDER','N_LENDERS','JUMBO','MC','HELD','ALTERNATIVE'),with=F]
  
  
  demand.data <- market.data[,c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','market.log_income.mean','market.log_income.std',
                                'market.log_price.mean','market.log_price.std','market.conforming_loan_limit','market.income_price.corr',
                                'market.log_income.mean.BAR','market.log_price.mean.BAR','MARKET_SIZE','MARKET_SHARE_DATA','RATE','market.log_loan_size.mean','LENDER','JUMBO','j','RATE_INST','delta','xi'),with=F]
  
  merged.cf.data <- merge(supply.data,demand.data,by=c('MARKET_ID','YEAR','CBSA','PURPOSE','FICO','j','LENDER','JUMBO'))
  
  
  # If Jumbo securitization...
  if(JUMBO_SEC_RATE < 999) {
    nb.rows <- merged.cf.data[LENDER != 'b']
    nb.rows[,j := j * 100]
    nb.rows[,JUMBO := 'j']
    nb.rows[,ALTERNATIVE := 1]
    nb.rows[,MC := MC + JUMBO_SEC_RATE] # Idea is you can securitize under exactly the same process but at some spread over what you'd do if GSE.
    merged.cf.data <- rbind(merged.cf.data,nb.rows)
    merged.cf.data <- merged.cf.data[order(MARKET_ID,j)]
    
  }
  
  toReturn <- list(market.data = merged.cf.data,linear.model = linear.model,demand.pars = demand.pars,raw.draws = raw.draws)
  
  return(toReturn)
}


estimation = 'estimation_final'

CR = NULL
LIM = NULL
RGSE = NULL
NAME = 'baseline'

runCounterfactualsMain(estimation)
runCounterfactualsEquity(estimation)
runCounterfactualsJumboSec(estimation)

runCounterfactualsMain_by_Risk(estimation,risk = 'LOW')
runCounterfactualsMain_by_Risk(estimation,risk = 'HIGH')